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Summary of Research 


1.0 Introduction 

Research in the FY01 focused on the analysis and development of enhanced algorithms 
for unsteady aerodynamics and chemically reacting flowfields. The research was performed in 
support of NASA Ames’ efforts to improve the capabilities of the in-house CFD code— OVER- 
FLOW. Specifically, the research was focused on the four areas given below. 

1. Investigation of Stagnation Region Effects. The preconditioning formulation, that is used for 
low speed aerodynamics computations in the OVERFLOW code, is known to be susceptible to 
instabilities in stagnation regions. In our research, we have employed perturbation analysis to 
shed light on the underlying causes of these instabilities. Specifically, it is evident that the insta- 
bilities are driven by the presence of large pressure fluctuations and it is necesary to modify the 
preconditioning formulation to account for such fluctuations. Based upon this insight, we have 
derived several modified preconditioning formulations and specific recommendations have been 
made for implementation in the OVERFLOW code. 

2. Unsteady Preconditioning Dual-Time Procedures. We have aided NASA researchers in 
implementing the unsteady preconditioning algorithm into the OVERFT.OW code. In addition, a 
multiple scales preconditioning formulation has been developed to address more complicated 
unsteady problems. This is also based upon using separate preconditioning formulations for the 
artificial dissipation terms and for the time-stepping. Detailed assessments of the enhanced algo- 
rithm need to be carried out and are currently underway. 
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3. Dissipation Formulation for Combustion. Various formulations for incorporating source 
terms into the artificial dissipation terms have been catalogued and evaluated. Appropriate forms 
for implementation into the OVERFLOW code have been identified. 

4. Time-Stepping Methods for Combustion. We have assisted NASA researchers in the imple- 
mentation of combustion capability in the OVERFLOW code. In addition, several methods of 
integrating the fluids equations with chemical source terms have also been evaluated. These 
include the fully implicit procedure, the Strang splitting method and the consistent splitting meth- 
ods. The fully implicit procedure has been found to be efficient and accurate but lacks robustness. 
On the other hand, the Strang splitting procedure is robust and reasonably efficient, but it is inac- 
curate near the combustion region. Finally, the consistent splitting method is robust and accurate, 
but suffers from inefficiencies. However, it is readily adaptable to the diagonalized solution proce- 
dure that is the solution algorithm of choice in the OVERFLOW code. Thus, the latter approach 
appears to be a promising avenue. Further studies are necessary to verify these findings. 

In the following sections, the specific research accomplishments in each of the above areas 
are discussed in greater detail. In addition, recommendations for further investigations are offered. 


2.0 Stagnation Region Effects 

Time-marching algorithms typically utilize preconditioning techniques to remain effective 
at low speeds or Mach numbers [1]. Although there are a wide variety of preconditioning tech- 
niques, all of them involve selective alteration or scaling of the physical time-derivatives in order 
to condition the system eigenvalues for all flow conditions. Specifically, all the methods essen- 
tially involve scaling of the pressure time-derivatives in the system by a term that is inversely pro- 
portional to the Mach-number squared. This scaling is crucial both for maintaining accuracy of 
the discrete formulation as well as for obtaining efficient convergence rates of the iterative solu- 
tion procedure. However, in some instances such as in the stagnation region, this scaling has the 
adverse effect of loss of robustness. Although this problem has been widely recognized, it is not 
very well understood. In this section, we analyze the underlying causes for the loss of robustness 
and then use this insight to devise amelioratory procedures. 
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The preconditioning method utilized in the OVERFLOW code is a variant of the form 
orginally proposed by Choi and Merkle [2] and then later modified by Weiss and Smith [3]. It is 
closely related to the procedures utilized by Turkel [4], but it is different from that suggested by 
Van Leer et al. [5], As mentioned above, in spite of their differences, all of these methods possess 
similar properties and all of them share a degradation in overall robustness for certain types of 
problems. To understand the reason for these difficulties, it is instructive to re-visit the derivation 
of the Choi-Merkle/Weiss-Smith formulation, which is based upon a perturbation analysis of the 
Euler equations. Here, we follow the procedure outlined by Venkateswaran and Merkle [1], 


2.1 Perturbation Analysis 

We start with the one-dimensional Euler equations for simplicity and note that a similar 
development is possible for the multi-dimensional Navier-Stokes equations as well. We write the 
equations in vector form as: 
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Here, we have written the equations using a primitive variables set. Although these are not neces- 
sarily the variables employed in the solution algorithm, we will see that the primitive variables set 
play a crucial role in the development of the preconditioning system. Note that the matrix T in 
Eqn. (1) represents the exact Jacobian of the variable transformation and the equations remain 
precisely equivalent to the standard physical system. 
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To analyze the behavior of the physical system under low Mach conditions, we employ 
perturbation expansions. We begin with the momentum equation, which we write in the noncon- 
servative form for simplicity, 


du du dp n 

+ P U— + 7T- = 0 




dx 


( 2 ) 


We next nondimensionalize the equations to facilitate order of magnitude comparisons 
between the terms. We introduce the following reference scales for the variables. 


L, u r , 



r 


Here, L represents a characteristic length scale, u r is the reference velocity and x r is the refer- 
ence time scale (to be defined later). p r is the reference pressure, while p,. is the reference mix- 
ture density. The nondimensionalized version of Eqn. 2 is then given by: 
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where the overbars indicate nondimensional quantities. 

We note that the coefficient of the pressure gradient term in the above equation becomes 
very large in the limit of low speed flows. To insure that the pressure gradient is always balanced 
by the convective terms, we define the small parameter, 

2 
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and consider the limiting form of Eqn. 3 as this parameter goes to zero. We then expand the pres- 
sure in a power series, 

P = Pq + ZP\ + ••• ( 5) 


and substitute it into Eqn. 3 to get 
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A more complete procedure would use analogous expansions for all the variables, but the results 
show that only the zeroth order quantities of the remaining variables appear in the final equations. 
Consequently, to minimize the algebra, we perturb only the pressure. 

The purpose of the time-derivatives in a time-marching scheme is to drive the equations to 
the desired steady state solution. To insure that this process is efficient, we select the characteristic 
time-scale such that the time-derivatives are of the same magnitude as the convective terms. Spe- 
cifically, in Eqn. 6, the coefficient of the time derivative must be of order unity so that 
1 = L/u r . This requirement clearly implies that the appropriate time scale is the convective 

time scale of the fluid particle. In that case, there remains no term to balance the 1 /e term in the 
pressure gradient term, leading us to conclude that, 


d _Po 

dx 


= 0 


(7) 


Moreover, for most problems, p q is fixed by the boundary condition and, hence, this quantity can 
be taken to be independent of time as well. In other words, p 0 is a constant and all changes in the 
pressure occur only in the first-order pressure. Note that this result makes an assumption on the 
size of the pressure fluctuations in the flowfield— i.e., the pressure fluctutions are of 0(e) . In fact, 
there may be situations wherein this stipulation may not hold, potentially leading to difficulties. 
We will discuss such exceptions later. 

Using Eqn. 7, the zeroth order momentum equation therefore takes the form, 


_ dU -_du d P 1 
9 dx + ?U dx + dx 


= 0 


( 8 ) 


Note that the £ term in the pressure gradient has been cancelled by the l/£ term that multiples 
it. All the terms in Eqn. 33 are clearly of order unity. In the time-marching framework, it is evi- 
dent that this equation provides an adequate means of updating the velocity u . However, the pres- 
ence of the first-order pressure p | implies that we must have a viable way of updating it from the 
continuity and energy equations. Accordingly, we consider these equations next. 
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To analyze the scaling of the continuity equation, we again non-dimensionalize and intro- 
duce the perturbation expansion, 
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The above equations has the time-derivative of the first-order pressure term present. However, this 
time -derivative become vanishingly small in the limit of e — > 0 . Thus, in this limit, there is no 
means of reliably updating the first-order pressure. In practice, the problem manifests itself in 
slow convergence of the time-marching procedure at low speeds. We further note that the situation 
is precisely the same as the singularity encountered in the incompressible limit. 

One way of alleviating the singularity problem is to alter the time-marching system such 
that it remains well-conditioned at all speeds. In the present context, this may be achieved through 
the introduction of a pseudo-property term in the pressure time -derivative such that it becomes 
order unity. Thus, 


ap_' 
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( 10 ) 


With this scaling, the time-dependent system is rendered suitable for time-marching computations 
at all flow speeds. This scaling procedure is generally referred to as preconditioning and, in the 
incompressible limit, it precisely corresponds to the artificial compressibility formulation [6]. 

As noted earlier, the preconditioning scaling of Eqn. 10 is contingent upon the assumption 
inherent in Eqn. 7. In problems involving pressure variations that are larger than 0(e) , one may 
therefore expect difficulties. By enforcing the requirement that p 0 is a constant, the first order 

pressure may become unnaturally large. In turn, the preconditioning scaling of the time-deriv- 
ative in the continuity equation may render this term large rather than of order unity (as it is 
designed to do). As a result, unstable solution behaviour may be encountered, leading to the 
observed lack of robustness. A common situation wherein this probem is encountered is in flows 
where the pressure level is not known a priori. For instance, in choked nozzle flows, the mass flow 
is often specified as an inflow boundary condition and the chamber pressure is determined as a 
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result of the solution. In such cases, large transients in the thermodynamic pressure {p 0 ) may 

occur, contributing to potential instability when the standard preconditioning scaling (in Eqn. 10) 
is adopted. Stagnation point flows are likewise associated with strong local pressure gradients, 
leading to pressure fluctuations, which may trigger local instabilities. 


2.2 Mitigating Procedures 

The perturbation analysis sheds light into the fundamental mechanics of the instabilities 
that arise due to large pressure fluctuations. Further, this insight suggests various avenues for mit- 
igating the problem. We discuss three potential methods below. 


2.2.1 Update Equation for the Zeroth Order Pressure 

The fundamental problem arises because of the assumption of a constant zeroth order 
pressure. This suggests that a possible solution would be to allow Pq to be time-varying. Of 
course, this means that an additional equation has to be introduced to provide a means of updating 
p Q , independent of /?, . A temporal equation for the mean pressure level in a chamber can be 

obtained by treating the region as a zero-dimensional volume (i.e., by using a lumped mass 
assumption). Using global conservation of mass, we get. 


3p d _Po 
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Augmenting the preconditioned system with the above equation for the thermodynamic pressuie 
allows the first order pressure to represent only the dynamic pressure. Thus, large pressure 
changes are handled by Eqn. 1 1 , while small pressure fluctuations are handled by the conventional 
preconditioned continuity equation, Eqn. 10. This formulation is very simple to implement; how- 
ever, its disadvantage is that it is a global “fix” and not a local one. Thus, it is well equipped to 
handle problems wherein the base pressure level is not known at the start of the computation, but 
it is not designed to deal with local pressure disturbances such as those encountered in stagnation 
regions. 
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2.2.2 Modified Preconditioning Scaling 

A second approach is to modify the scaling of the pressure time derivtaive in the continu- 
ity equation to account for the size of the local pressure fluctuation. In other words, the scaling 
suggested in Eqn. 10 is altered in such a manner that the pressure update is order unity for all 
magnitudes of Mach number and pressure fluctuation. Thus, 
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Here, when the pressure fluctuations are large, they are used to select the scaling of the time- 
derivative, until a point is reached, wherein the physical property derivative is itself well-condi- 
tioned. At that juncture, the system is no longer preconditioned. The advantage of this formulation 
is that it can be implemented on a local basis and is therefore be particularly effective for prob- 
lems involving stagnation regions. In fact, other researchers [3, 5] have successfully employed 
similar procedures for stagnation point flows. 

2.2.3 Linearized Preconditioning Procedures 

Yet another approach may be devised by making the observation that the instability arising 
from preconditioning scaling in the presence of large pressure fluctuations is fundamentally non- 
linear in nature. Therefore, a simple but effective way of handling the problem is to perform the 
preconditioning scaling after linearization rather than at the non-linear level. Thus, linear conver- 
gence is ensured even if the local pressure fluctuations are large. The linearized preconditioned 
scheme takes the following algorithmic form, 



Note that the formulation has been cast in a dual-time framework, wherein At is the non-linear 
time step and x is the linear time step. The preconditioning scaling is introduced only in the linear 
time-level, while the physical time-derivatives are retained in the non-linear level. Note that this 
formulation is completely general. In practice, it would be often possible to select the non-linear 
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time-step to be infinity, thereby leading to a single-time-step scheme. More details on the linear- 
ized preconditioning scheme is given in Ref. [1]. These studies and additional results will also be 
summarized in a forthcoming article to be presented at the AIAA Fluid Dynamics Conference [7], 


3.0 Unsteady Preconditioning Formulation 

Low Mach number and other stiff problems require the use of preconditioning scaling, 
which alters the form of the time-derivatives used in the time-marching algorithm. This means 
that the system is no longer time accurate. Unsteady computations therefore need to be performed 
within a dual time-stepping framework, similar in form to that presented in the previous section 
(Eqn. 13). The equations are marched in physical time and, at each physical time-step, the 
pseudo-time derivative is used to iterate the solution to convergence. The pseudo-time derivative 
may be altered or preconditioned to maximize convergence and accuracy. We point out that the 
preconditioning may be introduced at either the non-linear or the linear level, but we use the 
former procedure because it is customary to do so [1,8]. 

The OVERFLOW code has dual-time iterative capability for unsteady computations as 
well as preconditioning capability for steady problems. However, the original version of the code 
did not possess the ability to apply preconditioning for unsteady computations. The major task of 
this research was to provide NASA researchers with guidance regarding the implementation of a 
preconditioned dual-time procedure for unsteady computations. The main challenge in this regard 
is the tailoring of the dual time procedure to the diagonalized algorithm that is used in the OVER- 
FLOW code. We have followed the procedure of Merkle and Buelow [9] to accomplish this task. 
Secondly, the unsteady preconditioning procedure [1,8,9] was orginally devised for handling 
problems, characterized by a single unsteady time scale. However, this is inadequate for complex 
problems, wherein multiple local time scales may co-exist in a single problem. Accordingly, a 
second aspect of the research is to adapt the preconditioning strategy to handle multiple time 
scales. In the following sub-sections, we describe the two aspects of the research. 


3.1 Diagonalized Unsteady Preconditioning Formulation 

The preconditioned dual-time scheme may be expressed as: 
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Q . Q n + '-Q" . dE n+ ' , dF n+l = 
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where Y is the preconditioning matrix. This matrix premultiplies the pseudo-time derivative, 
while the other time-derivative retains its physical form. At each physical time-level, the precon- 
ditioned time-iterative procedure is exercised until a specified “convergence” level is attained. At 
that point, the solution is advanced to the next physical time level and the iterative procedure starts 
over. The choice of the preconditioning matrix is therefore decided as per the requirements of 
convergence and accuracy of the pseudo-iterations. The importance of accuracy is due to the 
dependence of the form of the artificial dissipation terms on the preconditioning matrix and will 
be addressed in greater detail in the following sub-section. 


The preconditioned implicit ADI version of the dual-time scheme may be written as: 
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As mentioned earlier, the main solution algorithm in the OVERFLOW code is the diago- 
nalized ADI procedure. Equation 15 thus needs to be adapted to the diagonalized procedure. Fol- 
lowing Buelow and Merkle [9], we get the following algorithm: 
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Now, introducing the modal matrices, we may diagonalize the LHS operators: 


M (l + AlA x ^M-jM^l + AxA y ^y-\Q k+ '-Q k ) = -A tS"'/ (17) 


The extension to 3D is straightforward and is not given here. 

3.2 Multiple Scales Preconditioning Strategy 

The definition of the preconditioning matrix for unsteady computations depends not only 
on the flow Mach number but also the choice of physical time scales. The choice is reflected in the 
precise definition of the pseudo-property term in the preconditioning matrix, which takes the fol- 
lowing form [1]: 

^ 2 - = Min 

dp 
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1 (nAt) 
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ap 

’d P 


(18) 


It can be observed that for low-frequency problems, the standard choice of the inverse-velocity 
squared is employed. However, for high frequency problems, a new unsteady velocity scale is 
used in the definition of the preconditioning parameter. 

The unsteady velocity scale in Eqn. 18 is computed from a global length scale and the 
physical time step size and is therefore a global velocity scale. Thus, it is appropriate only for 
problems that are characterized by a single unsteady frequency. In complex physical problems, 
multiple unsteady scales may be present. In such cases, the above definition needs to be modified 
in some fashion. However, one must bear in mind that preconditioning influences both the effi- 
ciency and accuracy of the computation and, in unsteady computations, the appropriate choice of 
preconditioner for efficiency is often different from the appropriate choice for accuracy. 

To accomodate such conflicting demands, we propose a multi-scales preconditioning strat- 
egy that is based upon employing different sets of preconditioners for the time-derivative and the 
artificial dissipation terms. Thus, we have 
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where T, is the time-derivative preconditioner chosen for efficient convergence properties and 

r„ is the dissipation preconditioner chosen for accuracy of the flux formulation. The selection of 

the pseudo-properties in the two preconditioning matrices will depend upon the specific con- 
straints for achieving efficiency and accuracy. Accordingly, 
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Note that the time-derivative preconditioner has the same form as in Eqn. 1 8, while the dissipation 
preconditioner uses difrerent local scales to dertemine the unsteady velocity scale. Exactly how 
this scale should be defined on a local basis remains a research issue, although it would appear 
that local length and time scales could be prescribed a priori from a knowledge of overall physics 
of the flow. Additionally, it is possible that the scheme may be more readily implemented using an 
AUSM-style dissipation procedure as described in Ref. 10. 

4.0 Artificial Dissipation Formulation for Combustion Problems 

There is currently an effort to bring in multi-species combustion capability into the 
OVERFLOW code. We are closely working with NASA researchers in this endeavour. Further, 
we are focusing on the new challenges associated with source terms, particularly those arising 
from stiff chemistry. Overall, the research is focused on two areas — namely, artificial dissipation 
formulation and time-stepping schemes. In the present section, we address techniques to modify 
the artificial dissipation formulation to account for dominant source terms. Specificially, we cata- 
log various approaches that have been proposed by researchers in the field. Time-stepping issues 
are dealt with in the following section. 

The design of artificial dissipation terms in source-dominated computations may be con- 
trolled by two sets of circumstances, namely, well-resolved sources and under-resolved sources. 
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In the former instance, all the gradients induced by the presence of source terms are well-resolved 
by the computational grid. In this case, the artificial dissipation terms should be selected so as to 
provide the minimal amount of damping that is necessary to suppress odd-even oscillations in the 
flowfield. The second instance arises when the computational grid is too coarse to resolve the 
source-induced gradients in the flowfield. Such problems arise frequently in chemically reacting 
flows, wherein some of the chemical time scales are extremely small, leading to very rapid relax- 
ation of the solution (or part of the solution) to equilibrium. In such cases, it is typically impracti- 
cal to fully resolve the source terms and, therefore, the computational solution is characterized by 
a discontinuity that separates two equilibrium states. The formulation of the artificial dissipation 
function must then be to capture these discontinuities accurately, i.e., they must ensure the correct 
location of the discontinuity, the correct equilibrium states on either side of it and solution mono- 
tonicity in its vicinity. 

The situation is similar to shocks, that arise from the solution of the Euler equations of 
compressible flow. In the absence of source-terms, several well-known methods have been 
devised to capture such discontinuities. These include the artificial dissipation methods of Jame- 
son and Turkel [11], the flux-vector splitting methods of Steger and Warming and Van Leer, the 
approximate Riemann methods of Roe [12] as well as the more recent AUSM method of Liou 
[13], The adaptation of these methods to account for the presence of source terms has been con- 
sidered by several researchers in different contexts. This includes some early work by Roe [14], 
and the more recent contributions of Pember [ 1 5], Bermudez and Vasquez [16], Jin and Lever- 
more [17], Mohanraj et al. [18], Yu and Chang [19] and Sussman [20]. While some of these 
papers are concerned primarily with under-resolved flows, many consider the well-resolved situa- 
tion as well. At any rate, we believe that the two sets of circumstances are closley related and can- 
not be studied in isolation. Many of these proposals are related and, in general, they all point to 
the need for a distributed treatment of the source terms. However, differences remain in the flavor 
and form of the implementation, and there is no general consensus amongst the authors. We will 
therefore review several of these contributions in this report with the aim of determining the right 
direction to pursue for further development. 

4.1 Modifications to the Roe Scheme 

Roe has proposed at least three different variants of his flux-differencing procedure to 
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accommodate the presence of source terms in the system of equations [12,14]. All three modifica- 
tions are somewhat related, but vary in the precise details. 

In Ref. 12, Roe considers quasi- ID flow with area variation. This corresponds to the 
source term vector being defined as: 


H = 



0 


( 21 ) 


The characteristic form of the quasi- ID Euler equations now becomes: 


In the steady state, 
the form: 


dW A dW ~ 

+ A^— = M l H = S 
dt dx 


these characteristic equations yield solutions that satisfy scalar equations 


( 22 ) 

of 


dw k _ s k 
dx \ k 
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Roe’s approach is to define a new Riemann variable based upon the above equation such that: 


8vv t ' = 8w/, 


s k hx 


( 24 ) 


and to re-define the flux-differencing procedure in terms of this new variable. The interfacial flux 
is then defined as: 


£. + 1 = ^E(Q L ) + E(Q R ))- ] -^\X k \bw k ffl k ( 25 ) 

In other words, the above procedure modifies the wave-strengths to account for the presence of 
source terms in the characteristic equations. 

The second approach used by Roe has been summarized by Barth [21]. The scheme 
appears similarly motivated — that is, it is modified to admit locally stationary solutions to account 
for the presence of the source term— but it is expressed in a different fashion. The numerical flux 
at the interface is modified to include the source term in the following fashion: 
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( 26 ) 


E i = E(Q l ) - ^\a\a-'H{Q l , Q„) 

' + 2 

which can, in turn, be written as: 

E i = \{E[Q L ) + E[Q K ))-] i \A\(Q s -QO-^\AA-'H(Q L ,Q l< ) (27) 

1 2 

Although Eqn. 27 is closely related to Eqn. 25, it is much easier to implement in practice. Also, 
note that the average source term may be computed by a trapezoidal rule integration. Further, 
Barth [21] makes the point that if the source discretization is carefully chosen, the overall scheme 
may be rendered higher-order accurate. 

For a first-order scalar system, it is instructive to examine Eqn. 25 further. Using simple 
averaging of the source term flux, we can write the full first-order discrete scheme as: 

|* + filial. A , + ^l^ + Tr (28) 

dt 25jc 2 dx 2 2 dr 

which involves an additional dissipation term, characterized by the source term. 

A third approach is outlined by Roe and Arora [14], where the authors adopt a more 
involved characteristics integration strategy to define the inteface fluxes. The method appears to 
be similar to the earlier approaches in principle and is not discussed further here. 


4.2 Bermudez and Vazquez Source Term Discretization 

Bemudez and Vazquez [16] find in their analysis of source-dominated flows that point rep- 
resentations of the source term have more limited stability and positivity constraints. They present 
a formal derivation to obtain a numerical source function, that parallels the development of the 
numerical inviscid fluxes that arise in flux-difference schemes. We summarize their development 
here. 

The overall discrete scheme is expressed as follows: 


E [ - E 

da . in i 

dt 8x 


1 

2 



( 29 ) 


where the numerical fluxes E , = E ](Q L , Q R ) , while the numerical source function, 

(± 2 J± 2 
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Hi = H i (x i _ ,, x t , x- + j, Qi ~ , j Qi, Q i+ i ) • The numerical fluxes are defined as given before for 
the Roe scheme, while the authors provide the following prescription for the numerical source 
function: 


H i (x i _ l ,x i ,x i+] ,Q i _ i ,Q i ,Q i+] ) = H L (x i _ v x i ,Q i _ ] ,Q i ) + H R (x i ,x i+l ,Q i ,Q i+l ) (30) 
where 

H L (x i _ l ,x i ,Q i _ l ,Q i ) = l -[I + \A\A~']H(x i _ l ,x l ,Q l _ l ,Q l ) 

H R (x i ,x i+i ,Q i ,Q i+ |) = ^[1 + \A\A~ ] ]H(x i ,x i+] ,Q i , Q i+{ ) 

Note the close similarity to Eqn. 26 of Roe’s second method. Bermudez and Vazquez recommend 
one of the following expressions for the average source values: 


H(x i _ ] ,x i ,Q i _ u Q i ) = x\H(x i _ u Q i _ ] ) + H(x i ,Q i )] 


(31) 


_ (x i _ ] + Xj Qi_]+ Qj 

//(*,_ „ *„ Q,_ ,, Q,) = H (-!— L -!-X 


(32) 


Substituting these expressions into Eqn. 30, we would obtain an averaged representation of the 
source terms. 


4.3 Other Approaches 

Jin and Levermore [17] offer two approaches to handle source term effects. Both methods 
are closely related to Roe’s approaches in that they seek to include the source term in the defini- 
tion of the interface numerical flux. The difference between the two approaches, in fact, lies in 
where the Taylor expansion is centered. In the piecewise steady method, similar to Roe, the Taylor 
expansion is performed about the nodal values, while in their so-called modified upwind scheme, 
the Taylor expansion is centered about the interface point. This subtle difference leads to philo- 
sophically related but distinct expressions for the source fluxes. Jin and Levermore prove that both 
their methods preserve proper solution behavior in the asymptotic limits of small, intermediate 
and large source terms. 

Mohanraj et al. [ 1 8] adopt a two-pronged strategy in addressing source term effects. They 
adopt Roe’s first strategy (Eqn. 25) and augment it by further using a distributed source represen- 
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tation similar to that suggested by Bermudez and Vasquez. They report improvements for several 
source-dominated flows. Further, they catalog higher-order source term discretizations as well in 
their paper. 

Finally, Yu and Chang [19] use the so-called space-time methodology to construct fluxes 
in the presence of source terms. They argue that unlike conventional flux-difference methods, the 
space-time method naturally introduces a distributed evaluation of the source term, which yields 
the correct physical behavior for source-dominated flows. 


5.0 Time-Stepping Schemes for Combustion Problems 

Over the years, several methods have been employe by researchers to compute reacting 
flow fields. Overall, we may classify the methods as being direct methods or splitting methods. In 
direct methods, the source term is treated fully implicitly, while in splitting methods, typically, the 
source term is treated using an ODE solver, while the convection-diffusion terms are treated in a 
conventional manner using a PDE solver. The most popular splitting method is the one due to 
Strang, but many variants are possible, among which we consider one that we refer to as the con- 
sistent splitting method. 


5.1 Direct or Fully Coupled Scheme 

In direct or fully coupled methods, the chemical source terms and the transport terms are 
solved simultaneously. In one dimension, the direct scheme may be written as: 


I + Ax D + Ax^ 

dx 


ias = -Axfff 


(33) 


where the residual R contains the convection diffusion terms, while the source term vector H 
contains the chemical source terms. 

Direct methods are often considered the best method because of its unconditional stability. 
In fact, this statement is not strictly true as discussed in detail in Ref. 22. Direct or implicit treat- 
ment of negative sources (or sinks) is unconditionally stable; however, implicit treatment of posi- 
tive sources is only conditionally stable. In combustion computations, it is normal to encounter 
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positive source terms during the initial non-linear stages of the computation, while, in the later lin- 
ear stages, only sink terms are encountered. Thus, full implicit treatment may lead to instabilities 
initially, but is asymptotically stable in the linear convergence process. 

A further difficulty of direct schemes is that they are not ideally suited for implementation 
in diagonalized implicit schemes. This is due to the fact that the source term Jacobian cannot be 
readily diagonalized. Thus, it is customary to solve the species equations (that contain the source 
terms) uncoupled from (or, more precisely, loosely coupled to) the fluids equations, which in turn 
may lead to additional loss of robustness and efficiency. For these reasons, we will consider the 
properties of the Strang splitting and consistent splitting schemes next. 


5.2 Strang Splitting Scheme 

In splitting methods, the integration of the equations are carried out in two steps, namely a 
chemistry integration, which requires the use of an ODE solver, followed by a transport integra- 
tion, which utilizes a conventional PDE solver. In fact, since the transport step does not contain 
the source term Jacobian, it is well suited for use with a diagonalized implicit solver. We may 
express the Strang scheme in the following fashion: 


Q 



Ax 


Ax 


J H r (x)dx 


At 


Q n + 1 - Q* + dE n+{ 

Ax dx 


= 0 


( 34 ) 


( 35 ) 


where we have first carried out a chemistry integration step followed by a transport-step. 

The Strang scheme is robust because the ODE solver automatically adapts the integration 
time-step according to the stability constraints of the problem. In other words, for positive 
sources, the ODE solver reduces the time-step size to insure stability and gradually increases it as 
the sink terms become dominant. Further, the procedure is efficient because it allows the use of 
the diagonalized algorithm for the transport step. However, the Strang splitting scheme does have 


Algorithmic Enhancements for Unsteady Aerodynamics and Combustion ApplicationsDecember 29, 2001 


18 



the disadvantage of compromising the accuracy of the computation. This is readily seen when the 
scheme is written in delta form in the “sink” limit: 


[I + AtD] 


I A 

/ + Ax^ 

Ox. 


A G = (36, 


It is evident that, upon convergence, the RHS operator retains a splitting error in the final solu- 
tion. In fact, this error is proportional to the source term Jacobian. Because of the relatively large 
magnitudes of the source term eigenvalues, this error term can be large, leading to serious inaccu- 
racies in the converged result. Example computations showing the effect of this inaccuracy are 
given in Ref. 23. 


5.3 Consistent Splitting Scheme 

It is desirable to combine the robustness of the splitting method with the accuracy of the 
direct method. This is accomplished in the following manner: 


Q 



Ax 



(37) 


Q n + _ 1 - Q_ + dE n + 1 
Ax dx 



(38) 


Note all terms in the governing equation are present in both chemistry and transport steps. In the 
chemistry step, the source term is advanced, while the transport terms are lagged. In the transport 
step, the transport terms are advanced, while the source term is evaluated at the intermediate * 
level. For this reason, this scheme is referred to as the consistent splitting method. 

The general accuracy of the formulation is readily verified by writing the scheme in the 
delta form: 


[/ + AxD] 


, A 

I + Ax^ 

Ox 


. fdE n 

Ae = - A T Tx ~ H ' 


( 39 ) 
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Note that there are no splitting errors introduced into the residual and, therefore, upon conver- 
gence, the scheme will converge to the proper solution. Further the consistent splitting scheme 
preserves the robustness property in the same manner as the Strang splitting scheme. However, 
there is a new source of concern, viz., the splitting error on the LHS. This error is akin to the 
approximate factorization error in ADI schemes and introduces a time-step limitation. By itself, 
the splitting error does not materially impact the choice of the time-step. However, when the con- 
sistent splitting scheme is combined with the ADI (or other factorization) method for the transport 
operator, the combined effect of the splitting error can introduce a significant limitation of the 
time-step size. 

The optimal time-step selection of the ADI scheme is well understood [24], The presence 
of the source term operator however introduces further complexity. Detailed investigation of these 
effects are summarized in Ref. 23 and are not given here. The final optimal time-step selection is 
given as follows: 

NDT = 2ndLargest(CFL CFL , CFL. STN ) (40) 

a V 4 

where NDT stands for non-dimensional time-step size and is usually around 10, CFL X stand 

for the CFL numbers and have their usual definitions, and STN is the source term number, 
defined using the largest eigenvalue of the source Jacobian. Thus, the optimal time-step should be 
based upon the second largest non-dimensional time-step in the problem, rather than the maxi- 
mum or the minimum. While this result is general, its implementation in a coupled system is not 
always straightforward and further studies are necessary to refine it. 


6.0 Summary 

Research in FY 2001 has touched on a wide array of topics. The major goal has been to 
assist NASA researchers in the implementation of unsteady preconditioning and combustion 
capability in the OVERFLOW code. Further, research has been pursued in areas related to the 
above implementation. Specifically, the research falls into four categories: 
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1. Stagnation Region Effects. Perturbation analysis of the low Mach number limit of the Euler 
equations has been carried out to shed light on the underlying causes of difficulties. In particular, 
the analysis has revealed that large pressure fluctuations can lead to unstable behavior in the pres- 
ence of preconditioning. It is believed that this is the fundamental reason for the cause of instabil- 
ities near stagnation regions. The perturbation analysis has also pointed to several mitigating 
procedures, namely, updating the zeroth order pressure independently, using a modified precondi- 
tioning choice and applying preconditioning at the linear level. Of these, the modified precondi- 
tioning choice which considers the magnitudes of the local pressure fluctuations appears to be the 
most promising procedure. Implementation and testing in the OVERLFOW code are in progress. 

2. Unsteady Preconditioning Formulation. A diagonalized preconditioning formulation has 
been implemented in the OVERFLOW code by NASA researchers and is currently being evalu- 
ated for a variety of aerodynamic problems. In addition, we have developed a multiple scale pre- 
conditioning formulation, which is a refinement of the present method. Following further testing 
and verification, the improved method will be incorporated into the OVERFLOW code. 

3. Dissipation Formulation for Combustion. Methods of modifying the artificial dissipation to 
account for dominant chemical source terms have been reviewed. Overall, the different 
approaches are related. They essentially involve two aspects, namely the modification of the inter- 
face numerical flux formulas and the representation of an averaged source term. Further, the deri- 
vation of higher order representations of the system including the source term needs to be studied 
further. The most promising approaches will then be compared and evaluated for implementation 
in OVERFLOW. 

4. Time-Stepping Schemes for Combustion. Three different time-integration techniques for 
combustion problems have been compared and contrasted. The direct (or fully implicit) scheme 
has been observed to be efficient and accurate but lacks robustness. The Strang splitting scheme is 
robust and efficient, but is inaccurate under some circumstances. Finally, the consistent splitting 
scheme is robust and accurate, but is inefficient under some circumstances. Nevertheless, it 
appears that the consistent splitting scheme is well-suited for implementation within a diagonal- 
ized framework. Therefore, it shows promise for implementation in OVERFLOW. Moreover, it 
may be possible to combine the method with schemes that are already used in OVERFLOW. 
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